#ifndef __CRAZYFISH_MATRIX__
#define __CRAZYFISH_MATRIX__

#include <iostream>
#include <new>
#include <vector>
#include <exception>
#include <cmath>
#include <limits>

/// Set the machine epsilon.
const double eps = std::numeric_limits<double>::epsilon();

template <typename T>
class SquareMatrix;

template <typename T>
std::vector<T> operator/ (const std::vector<T>& lhs, const SquareMatrix<T>& rhs);

template <typename T>
class SquareMatrix
{
private:
    T* matrix_elements;
    int dimension;
    void allocMem(int _d);
    void releaseMem();
    void exitAbnormal(int _err_code) const;
    void copy_data(const T* _f, T* _t, int _d);
    
public:
    SquareMatrix()
	{
	    matrix_elements = nullptr;
	}
    SquareMatrix(int _n);
    SquareMatrix(const SquareMatrix<T>& rhs);
    ~SquareMatrix();
    SquareMatrix& operator=(const SquareMatrix<T>& rhs);
    T& operator()(int _i, int _j);
    std::vector<T> invMul(const std::vector<T>& rhs) const;
    
    friend
    std::vector<T> operator / <> (const std::vector<T>& lhs, const SquareMatrix<T>& rhs);
};

void exitAbnormal(int error_code);

void exitAbnormal(int error_code)
{
    switch (error_code)
    {
    case 1:
	std::cerr << "The system is linear dependence." << std::endl;
	break;
    case 2:
	std::cerr << "The dimensions of the vector and the matrix are not match." << std::endl;
	break;
    default:
	break;
    }
    std::terminate();
};

template <typename T>
void SquareMatrix<T>::copy_data(const T* _f, T* _t, int _d)
{
    for (int i = 0; i < _d; i++)
	_t[i] = _f[i];
};

template <typename T>
void SquareMatrix<T>::exitAbnormal(int _err_code) const
{
    switch (_err_code)
    {
    case 1:
	std::cerr << "index out of range." << std::endl;
	break;
    case 9:
	std::cerr << "no enough memory." << std::endl;
	break;
    default:
	break;
    }
};

template <typename T>
void SquareMatrix<T>::releaseMem()
{
    if (matrix_elements != nullptr)
	delete [] matrix_elements;
    matrix_elements = nullptr;
};

template <typename T>
void SquareMatrix<T>::allocMem(int _d)
{
    try
    {
	matrix_elements = new T [_d];
    }
    catch (std::bad_alloc& no_memory)
    {
	std::cerr << no_memory.what() << std::endl;
	exitAbnormal(9);
    }
};

template <typename T>
T& SquareMatrix<T>::operator()(int _i, int _j)
{
    if (_i < 0 || _i > dimension - 1 || _j < 0 || _j > dimension - 1)
	exitAbnormal(1);
    return *(matrix_elements + _i * dimension + _j);
};

template <typename T>
SquareMatrix<T>& SquareMatrix<T>::operator=(const SquareMatrix<T>& rhs)
{
    int n = rhs.dimension;
    if (matrix_elements != nullptr)
	releaseMem();
    allocMem(n * n);
    copy_data(rhs.matrix_elements, matrix_elements, n * n);
    dimension = n;
    return *this;     
};

template <typename T>
SquareMatrix<T>::SquareMatrix(const SquareMatrix<T>& rhs)
{
    int n = rhs.dimension;
    allocMem(n * n);
    copy_data(rhs.matrix_elements, matrix_elements, n * n);
    dimension = n;
};

template <typename T>
SquareMatrix<T>::SquareMatrix(int _n)
{
    allocMem(_n * _n);
    dimension = _n;
};

template <typename T>
SquareMatrix<T>::~SquareMatrix()
{
    releaseMem();
};

template <typename T>
std::vector<T> SquareMatrix<T>::invMul(const std::vector<T>& rhs) const
{
    int n = dimension;
    std::vector<T> x(n);
    SquareMatrix<T> A(*this);
    std::vector<T> b(rhs);
    if (n != rhs.size())
	exitAbnormal(2);
    /// A gaussian elimination procedure without pivoting.
    for (int k = 0; k < n - 1; k++)
    {
	/// Check the linear dependency of the system.
	if (std::fabs(A(k, k)) < eps)
	    exitAbnormal(1);
	for (int i = k + 1; i < n; i++)
	{
	    for (int j = k + 1; j < n; j++)
		A(i, j) -= A(k, j) * A(i, k) / A(k, k);
	    b[i] -= b[k] * A(i, k) / A(k, k);
	}
    }
    if (std::fabs(A(n - 1, n - 1)) < eps)
	exitAbnormal(1);

    /// Backward substitution.
    for (int i = n - 1; i >= 0; i--)
    {
	T sum = 0.0;
	for (int j = i + 1; j < n; j++)
	    sum += A(i, j) * x[j];
	x[i] = (b[i] - sum) / A(i, i);
    }
    return x;
};

template <typename T>
std::vector<T> operator/(const std::vector<T>& lhs, const SquareMatrix<T>& rhs)
{
    int n = rhs.dimension;
    std::vector<T> x(n);
    x = rhs.invMul(lhs);
    return x;
};

#else
/// DO NOTHING.
#endif
